% PERSPECTIVE TRANSFORMATION
% ==========================
%
% This source file illustrates the perspective transformation combined with
% image convolution. The outlined functionality can be achieved by using
% matlab's inbuilt functions, but that is not helpful for implementation
% in other programming environment. Therefore this example shows step-by-step
% approach using matlab's basic functions. 
%
% First step is the calculation of the transformation matrix:
%
% Define the corner points around the area you want to transform into rectangle
% The corners of the trapezoid used in the example are:
%
%     [-100,1]      [1100,1]
%        x              x
%
%
%          x         x
%       [1,1450] [1000,1450] 
%
% The transformation matrix is a classic 3x3 matrix known from all other point
% operations. The matrix has 9 elements. The 4 points are used in the equations 
% for calculating the x,y coordinates from u,v plane. That delivers 8 equations. 
% The ninth element in the transformation matrix - a33 is considered 1, while 
% that parameter is responsible for overall scale only.
%
% The equations are determined from the classic transformations:
%
%                          | a11 a12 a13 |
%     [x',y',w'] = [u,v,w] | a21 a22 a23 |
%                          | a31 a32 a33 |
%
% where the relationship x = x'/w' and y = y'/w'.
%
% The corresponding calculation:
%
%         x'     a11u + a21v + a31
%    x = ---- = -------------------
%         w'     a13u + a23v + a33
%
%         y'     a12u + a22v + a32
%    y = ---- = -------------------
%         w'     a13u + a23v + a33
%
% unlike the affine transformation the variable w' is not constant and may take different values
% for different u,v combinations. Therefore we need to consider it in the calculations.
%
% The equations can be expressed as:
%
%    x = a11u + a21v + a31 - a13ux - a23vx
%    y = a12u + a22v + a32 - a13uy - a23vy
%
% This equations are taken for 4 corner points of the enveloping trapezoid and the system 8x8
% emerges:
%
%    | u0  v0  1  0  0  0 -u0x0 -v0x0 |   | a11 |   | x0 |
%    | u1  v1  1  0  0  0 -u1x1 -v1x1 |   | a21 |   | x1 |
%    | u2  v2  1  0  0  0 -u2x2 -v2x2 |   | a31 |   | x2 |
%    | u3  v3  1  0  0  0 -u3x3 -v3x3 |   | a12 |   | x3 |
%    |  0   0  0 u0 v0  1 -u0y0 -v0y0 | . | a22 | = | y0 |
%    |  0   0  0 u1 v1  1 -u1y1 -v1y1 |   | a32 |   | y1 |
%    |  0   0  0 u2 v2  1 -u2y2 -v2y2 |   | a13 |   | y2 |
%    |  0   0  0 u3 v3  1 -u3y3 -v3y3 |   | a23 |   | y3 |
%
% which is a classic example for matrix based solution. Use premultiplication with inverse
% matrix A, which yields identity matrix on left side. By calculating the right side we get
% the desired result:
%
%    M.a=x
%    inv(M).M.a=inv(M).x
%    a=inv(M).x
%
% The result will deliver the coefficients a in a form of column vector, where the elements are
% arranged like this:
%
%  | a11 |
%  | a21 |
%  | a31 |
%  | a12 |
%  | a22 |
%  | a32 |
%  | a13 |
%  | a23 |
%
% ,which we rearrange into the transformation matrix:
%
%                          | a11 a12 a13 |
%     [x',y',w'] = [u,v,w] | a21 a22 a23 |
%                          | a31 a32 a33 |
%
%   NOTE: this example uses 'postmultiplication' syntax. It is equivalent to 'premultiplication'
%   syntax, which is frequently found in literature. The transformation matrix must be
%   transposed in order to swap for premultiplication syntax. 
%
% The test matrix based on the trapezoid described earlier:
A=[
	   1,   1,   1,   0,   0,   0,    100,      100;
	1000,   1,   1,   0,   0,   0, -1.1e6,    -1100;
	1000,1450,   1,   0,   0,   0,   -1e6,  -1.45e6;
	   1,1450,   1,   0,   0,   0,     -1,    -1450;
	   0,   0,   0,   1,   1,   1,     -1,       -1;
	   0,   0,   0,1000,   1,   1,  -1000,       -1;
	   0,   0,   0,1000,1450,   1,-1.45e6,-2.1025e6;
	   0,   0,   0,   1,1450,   1,  -1450,-2.1025e6];

% Trapezoid around the image
x=[-100;1100;1000;1;1;1;1450;1450];

% ***** Small Sample Imaging *****
% This is matrix for small sample image (testing the convolution):
A=[
	   1,   1,   1,   0,   0,   0,     22,       22;
    512,   1,   1,   0,   0,   0,-273408,  	 -534;
    512, 629,   1,   0,   0,   0,-262144,  -322048;
      1, 629,   1,   0,   0,   0,     -1,     -629;
      0,   0,   0,   1,   1,   1,     -1,       -1;
      0,   0,   0, 512,   1,   1,   -512,       -1;
      0,   0,   0, 512, 629,   1,-322048,  -395641;
      0,   0,   0,   1, 629,   1,   -629,  -395641];		
		
% Trapezoid around the small sample image
x=[-22;534;512;1;1;1;629;629];
xshow=[-22,1;534,1;512,629;1,629]; % Just for plotting purpose, not used in calculation

% Calculate the coefficients by solving the 8 linear equations.
iA=inv(A);
a=iA*x;

disp('Coeficients for forward mapping');
disp(a);

% Test the matrix for forward mapping:
% u=[1,1];
% for i=1:20:512
% 	u(1,1)=i;
% 	x=(a(1,1)*u(1,1)+a(2,1)*u(1,2)+a(3,1))/(a(7,1)*u(1,1)+a(8,1)*u(1,2)+1);
% 	y=(a(4,1)*u(1,1)+a(5,1)*u(1,2)+a(6,1))/(a(7,1)*u(1,1)+a(8,1)*u(1,2)+1);
% 	out(i,1)=i;
% 	out(i,2)=x;
% 	out(i,3)=y;
% endfor
% 
% % Show the resulting mapping ...
% disp("Forward mapping");
% disp(out);

% The outlined process is valid for so called forward mapping. That is a very useful 
% strategy however sometimes may be more convenient to use inverse mapping. Forward
% mapping traslates the integer inout into real number output. The inverse mapping 
% uses integer output matrix and maps it into real number in the input. That is easier
% to calculate in this application.
%
% Inverse mapping is the same process but the matrix of a-coefficients is inverted. 
% Create 3x3 matrix, calculate inverse matrix and set the new coefficients in the 
% a-vector again.
fwdA=[1,1,1;1,1,1;1,1,1];
fwdA(:,1)=a(1:3);
fwdA(:,2)=a(4:6);
fwdA(1:2,3)=a(7:8);

ia=a;
invA=inv(fwdA);
ia(1:3)=invA(:,1);
ia(4:6)=invA(:,2);
ia(7:8)=invA(1:2,3);

disp('Coeficients for inverse mapping');
disp(ia);

% Now we have a-vector based on inverse matrix.

% Test the matrix for inverse mapping:
% index=1;
% u=[-100,1];
% u=[-22,1];
% for i=-22:20:534
% 	u(1,1)=i;
% 	x=(ia(1,1)*u(1,1)+ia(2,1)*u(1,2)+ia(3,1))/(ia(7,1)*u(1,1)+ia(8,1)*u(1,2)+1);
% 	y=(ia(4,1)*u(1,1)+ia(5,1)*u(1,2)+ia(6,1))/(ia(7,1)*u(1,1)+ia(8,1)*u(1,2)+1);
% 	iout(index,1)=i;
% 	iout(index,2)=x;
% 	iout(index,3)=y;
% 	index++;
% endfor

% Show the resulting mapping ...
% disp("Inverse mapping");
% disp(iout);

%
%
% Determine the optimal interpolation function for signal reconstruction
% =======================================================================
%
% Filter: Convolution between a non-zero width pulse with the triangle kernel.
% -----------------------------------------------------------------------------
%
% This is a silly approach that considers the signal as a sequence of adjacent
% pulses. One can interpret this is double interpolation. First, nearest 
% neighbor interpolation is used, which yields the continuous, zip-zapped 
% signal. The width of a pulse is the distance between the pixels. This 
% approach establishes continuity in the sampled signal, which is important 
% for next assumption. The levels are  constant between pixels, and we know the
% intensity for arbitrary spatial index. Graphically can be expressed like
% this:
%
%             .-x-. 
%             |   |
% -x-.    .-x-'   '-x-.
%    |    |           |
%    .-x-.            '-x-.-x-
% 
% where the sampled values are actually denoted by 'x' symbol. One can split
% this signal into train of non-overlapping pulses. Each pulse can be convolved
% with various kernels. Triangle kernel and windowed sinc function kernel were
% tested in this demo.
% 
% This example shows that this rectangular pulse of non-zero width can be
% convolved with triangle function. Convolution yields non-linear curve. The
% further code shows how can be used the resulting curve for interpolation. 
%
% Helper functions and kernel functions
% --------------------------------------
%
% Test for integration - this function is difficult to integrate
function y = foo(x)
	y = x .* sin (1./x) .* sqrt (abs (1 - x));
endfunction

% Triangel kernel -1 - 1. It is described by two linear functions
function y=kernel_P(x)
	y=x.*-1.0 + 1.0;
endfunction

function y=kernel_L(x)
	y=x.*1.0 + 1.0;
endfunction

function y=kernel_2P(x)
	y=x.*-0.74 + 0.525;
endfunction

function y=kernel_2L(x)
	y=x.*0.74 + 0.525;
endfunction

% Calclucate manually the convolution of ideal pulse and triangel function.
idx=linspace(-0.5,0.5,128);
for i=1:128
	bottom=idx(i);
	top=idx(i)+1;
	if( top > 1.0 )
		top = 1.0;
	endif

	if( bottom < 0.0 )
		LEFT=quadcc("kernel_2L",bottom,0.0);
		RIGHT=quadcc("kernel_2P",0.0,top);
		I=LEFT+RIGHT;
	else
		I=quadcc("kernel_2P",bottom,top);
	endif

	ptemp(i)=I;
endfor

itemp=fliplr(ptemp);

temp=itemp;
temp(1,129:256)=ptemp;

trimmer=linspace(-1,1,256);

% Normalize the curve
if(min(temp) < 0.0)
	linconv=temp.+abs(min(temp));
else
	linconv=temp;
endif
linconv=linconv./(max(linconv)+min(linconv));

% Apply the ideal pulse filter to remove the aliasing artefacts.
% Trimming makes no sense, while the signal is not superposed, just stacked.
%trimmer=linspace(-1,1,256);
%for i=1:256
%	if((trimmer(i) < -0.5 ) || (trimmer(i) > 0.5))
%		linconv(i)=0;
%	endif
%endfor
%linconv=linconv./(max(linconv));

figure(1);
subplot(1,3,1)
plot(trimmer,temp),grid,title("Convolution with triangle");
subplot(1,3,2)
plot(trimmer,linconv),grid,title("Normalised Convolution");

% Test the convolution ...
pixel_1_intensity=100;
pixel_2_intensity=120;

for t=1:99
	position_after_transposition=(t/100);
	rPix=uint8(position_after_transposition * 128);
	lPix=uint8(rPix+128);
	reconstructed_intensity(t)=pixel_1_intensity*linconv(lPix)+pixel_2_intensity*linconv(rPix);
endfor

subplot(1,3,3);
plot(reconstructed_intensity),grid,title("Reconstructed intensity, interpolation is convolution between triangle and square pulse");

% end of silly experiment with convolution between triangle and non-zero-width pulse.
% =============================================================================================
%
% Filter: Convolute the samples with triangle function
% -------------------------------------------------------
% The samples are represented as zero-width pulses. 
%
%   x     x  x
%   |  x  |  |
%   |  |  |  |
% --'--'--'--'--
%
% The convolution is actually equal to convolution kernel. Therefore the
% convolution with triangle filter will produce linear interpolation between
% single pixels (samples) in the digitalised input.
% 
% Demo the linear regression:
for t=1:99
	position_after_transposition=(t/100);
	reconstructed_intensity_lin(t)=pixel_1_intensity+(((position_after_transposition-0)/(1-0))*(pixel_2_intensity-pixel_1_intensity));
endfor

% Interpolation with function generated by convolving windowed sinc with square pulse
% ------------------------------------------------------------------------------------
 
% Windowed sinc(x) function. The window is defined for the interval <-2pi,+2pi>,
% The windowing is important to scrap off the side lobes. This function is
% refitted for interval <-1,+1>. The convolution with the square will expand the
% interpolation function into interval <-2,+2>. The interpolation takes the
% central part of the convolution result from the interval <-1,+1> 

function y=sinc(x)
	y=sin(x)./x;
endfunction

sinc_filter=sinc(linspace(-2*pi,2*pi,256));
sinc_base=linspace(-2*pi,2*pi,256);
square=[zeros(1,64) ones(1,128) zeros(1,64)];

kernel_sinc=conv(square,sinc_filter); % convolution expands the size of the base interval to 511 elements
kernel_sinc=kernel_sinc./(max(kernel_sinc)-abs(min(kernel_sinc))); % Normalize the kernel so that intensities are preserved.

% Display the result
figure(2);
subplot(1,3,1);
plot(sinc_base,sinc_filter),grid,title("Windowed sinc function"),xlabel("Angle");

subplot(1,3,2);
sinc_base=linspace(-1,1,256);
kernel_base=linspace(-2,2,511);
plot(kernel_base,kernel_sinc,"1",sinc_base,square,"2",sinc_base,sinc_filter,"3"),
		grid,
		title("Convolution square signal with sinc"),
		xlabel("pixel grid");

for t=1:99
	position_after_transposition=(t/100);
	lPix=128+uint16(position_after_transposition * 128);
	rPix=256+uint16(position_after_transposition * 128);
	reconstructed_intensity_l(t)=pixel_1_intensity*kernel_sinc(rPix);
	reconstructed_intensity_r(t)=pixel_2_intensity*kernel_sinc(lPix);
	reconstructed_intensity(t)=pixel_1_intensity*kernel_sinc(rPix)+pixel_2_intensity*kernel_sinc(lPix);
endfor

subplot(1,3,3);
plot(reconstructed_intensity_l, "2", reconstructed_intensity_r, "3", reconstructed_intensity, "1"),
		grid,
		title("Interpolation with Windowed sinc function convolved with square pulse"),
		xlabel("distance between two pixels (percent)"),
		ylabel("intensity");

figure(3);
plot(reconstructed_intensity_lin, "1", reconstructed_intensity, "4" ),
		grid,
		title("Interpolation linear vs. sinc"),
		xlabel("distance between two pixels (percent)");

%
% PERSPECTIVE TRANSFORMATION & CONVOLUTION AGAINST VARIOUS KERNELS
% ----------------------------------------------------------------
% Read the test image, apply the transformation matrix, inverse perspective matrix for inverse
% mapping and use convolution to determine new intensity.
% The image is the same one as used for transformation implemented by native matlab's functions, 
% but minified to resoltuion 512x629 pixels. The corresponding trapezoid is given earlier in this
% sample.
% 

% Read the test image into matrix.
pkg load image
figure(4);
im=imread('sample.small.jpg');
imshow(im);

% minify the number of layers to one only - gra image.
g=rgb2gray(im);

% increase the image canvas avout 200 pixels. 100 left, 100 right.
rg=zeros(629,712);
rg(:,101:612)=g;
rg=uint8(rg);
 
% Show the grayscale image
figure(5);
imshow(rg);

% Perspective Transformation for x-dimension only. y-dimension is not relevant.
% Use the coefficients calculated in the first step, valid for inverse mapping. 
function x=ptrans(vc,u,v)
	ia=vc;
	x=(ia(1,1)*u+ia(2,1)*v+ia(3,1))/(ia(7,1)*u+ia(8,1)*v+1);
endfunction

% The convolution expands the spacial dimension 2x-1. The interpolation function
% is defined for the central part of the function. Convolution base was 256 stops,
% which results into 511 stops in the result. We extract 256 central stops, -128/+128
% from the cetral part.
function y=filter_sinc(kernel,value,intensity0,intensity1)
	pix_1_intensity=intensity0;
	pix_2_intensity=intensity1;
	position_after_transposition=value;
	lPix=128+uint16(position_after_transposition * 128);
	rPix=256+uint16(position_after_transposition * 128);
%	reconstructed_intensity_l(t)=pix_1_intensity*kernel_sinc(rPix);
%	reconstructed_intensity_r(t)=pix_2_intensity*kernel_sinc(lPix);
	y=pix_1_intensity*kernel(rPix)+pix_2_intensity*kernel(lPix);
endfunction

% ========== FORWARD MAPPING, PERSPECTIVE TRANSFORMATION =====================
% Process the input canvas pixel-by-pixel, use forward mapping. The canvas size of the output image
% is enlarged beforehand so that the input can be blown out into larger output. The forward transformation
% requires accumulator, usually much larger than the image. The accumulator must be present during
% the entire image is processed, theoretically we don't know where the output will exactly land. Input
% is coming in scanline order, but output is placed randomly in the output array, positioning is in
% real numbers. The accumulator must use a smart algorithm to accumulate the overlapping input.
erg=zeros(629,712);
erg=uint8(erg);
for r=1:rows(g)
	for c=1:(columns(g)-1);
		vx=ptrans(a,c,r);
		idx=vx+100;

		ratio=idx-floor(idx);

		% Accumulator. We have the whole image erg already and we check the existing value in it. 
		intensity=erg(r,floor(idx)); 

		intensity0=g(r,c); 
		intensity1=g(r,c+1);
		
		intensity0 = intensity + (intensity0 * (1-ratio));
		intensity1 = intensity1 * ratio;

		erg(r,floor(idx))=intensity0;
		erg(r,ceil(idx))=intensity1;
	endfor
endfor

figure(6);
imshow(erg);
hold on;
plot(xshow(:,1),xshow(:,2),"or");
imwrite(erg, 'sample.small.linear.size.exp.fwd.png', 'png');

% ========== INVERSE MAPPING, PERSPECTIVE TRANSFORMATION =====================
% Process the output canvas pixel-by-pixel. Use inverse mapping, the mate for the new image is increased
% for 100 pixes in both directions (100 left, 100 right). But we generally we need just single strip of the 
% output image. We walk it index by index (integers) and map in into the input image. We use linear regression
% (or more sophisticated convolution kernel) to get intensity value for a whole number output position.
% The inverse mapping is used for area defined inside the trapezoid only. The trapezoid elements are calculated 
% from forward transformation. The inverse mapping requires inverse matrix calculation from forward one.
% That operation is not really cheap, but the inverse mapping can then save memory and delivers better 
% results compared to forward, accumulator based mapping.
erg=zeros(629,712);
erg=uint8(erg);
for r=1:rows(rg)
%	disp(r);
	for c=1:columns(rg)
		idx=(c-101); % The integer-based output mate is defined for interval <-100, 612> px
		if(( idx > -22 ) && ( idx < 534 )) % The limitting trapezoid is <-22,534> px
			% nearest neighbor interpolation is simply created by rounding to
			% narest number - function uint16() in this case.
%			erg(r,c)=rg(r,(uint16(ptrans(ia,idx,r))+101)); 

			% Run transformation
			vx=ptrans(ia,idx,r);

			% Prepare data for approximation ...
			ratio=vx-floor(vx);
			intensity0=rg(r,(floor(vx)+101)); 
			intensity1=rg(r,(ceil(vx)+101)); 

			% linear regression - convolution with triangle kernel.
			erg(r,c)=intensity0+(ratio*(intensity1-intensity0));

			% Interpolation with "windowed sinc convolved with pulse signal".
			% The calculation is great deal inefficient compared to rounding,
			% and takes therefore way too much time. The interpolation is not
			% ideal, it produces funny moire effect, which needs to be filtered
			% after convolution again.
%			erg(r,c)=filter_sinc(kernel_sinc,ratio,intensity0,intensity1);
		endif
	endfor
endfor

figure(7);
imshow(erg);
hold on;
plot(xshow(:,1),xshow(:,2),"or");
imwrite(erg, 'sample.small.linear.size.exp.inv.png', 'png');
